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ABSTRACT 

We address the question of astronomical image processing from data obtained 
with array detectors. We define and analyze the cases of evenly, regularly, and 
irregularly sampled maps for idealized (i.e., infinite) and realistic (i.e., finite) 
detectors. We concentrate on the effect of interpolation on the maps, and the 
choice of the kernel used to accomplish this task. We show how the normalization 
intrinsic to the interpolation process must be carefully accounted for when dealing 
with irregularly sampled grids. We also analyze the effect of missing or dead 
pixels in the array, and their consequences for the Nyquist sampling criterion. 

Subject headings: methods: data analysis — techniques: image processing 

1. Introduction 

The creation of smooth two-dimensional maps from a series of samples measured at 
discrete points is a common problem in astronomical image processing. The goal is to create 
a smooth map with the best possible spatial resolution given a set of data sampled in two- 
dimensions. The solution is complicated by the fact that the data are often not sampled in 
a regular way even if the detector layout is regular. For example, telescopes may be scanned 
or dithered to map areas larger than the array, some instruments are unable to follow objects 
as they rotate on the sky, or an array itself may contain flaws (i.e., missing or dead pixels). 
The resulting two-dimensional sample pattern can often appear quite irregular. 
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Although the layout of most modem detector arrays (e.g., CCDs) can reasonably be 
approximated as generating evenly sampled grids (ESG) extending to infinity in all direc- 
tions, observations with these arrays will typically include a series of array translations and 
rotations. Additionally, one may want to combine multiple images of the same piece of sky 
in order to increase the signal-to-noise ratio. This requires that the images be registered so 
that the same area of the sky is being observed in each image. That is, any relative transla- 
tion and rotation of the array positions with respect to the sky must be taken into account 
when combining the images. Unless the translation and rotation operations are such that 
every pixel lies in a location previously occupied by another pixel then the resulting sample 
pattern is no longer an ESG. 

The processing of data from any case other than an ESG requires performing an inter- 



polati on. Some of these cases have been discussed by other authors (e.g.. lGranrath &: Lersch 



(119981 )). The interpolation (or smoothing) necessarily has an effect on the spatial resolution 
of the resulting map. For the case of any sampling pattern (ESG or otherwise) we wish to 
address the following two questions: 1) How does one choose an optimal kernel shape and 
size for the interpolation function? and 2) How does this kernel choice affect the spatial 
resolution of the resulting map? We will concentrate on the case of a Gaussian kernel. 

As mentioned above, the construction of maps from non-ESG sampled data is gen- 
erally done through an inte rpolation of the data using a smoothing kernel (see ^ and 



Lombardi fc Schneider! (120011 )). In this paper, we begin by reviewing the solution for ESGs 
using a technique based on Fourier transforms (§3]) and extend this technique in §l]to regu- 
larly spaced grids (RSGs), which are composed of relatively translated ESGs. In ^we use 
the tools developed for studying RSGs and ESGs to analyze irregularly sampled grids (ISGs) 
and explore the effects of missing samples in a map. 

Throughout this paper we present examples that reference the bolometer array used 
with SHARP. SHARP is a polarimeter module that is used in conjunctio n with the SHARC - 



II camera, which is deployed at the Caltech Submillimeter Observatory (iDowell et al.ll2003l ). 
For this, the 12 x 32 SHARC-II detector array is optically split into two 12 x 12 sub-arrays 
(and a section of 12 x 8 unused pixels), w hich image two orthogonal linear polarization 



components of radiation (INovak et al.ll2004l ). Although we concentrate on SHARP maps. 



the results are applicable to any detector array or sampling pattern. 
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2. Mathematical Definitions 

Before embarking on the analysis of the evenly sampled grid (ESG) we first introduce 
a set of definitions and functions that are central to the development of the subsequent 
sections. Given a function g (r), which is dependent on position r = xe^ + ye.yAe.x and Gy 
are the usual Cartesian unit basis vectors), we define the Fourier transform paiio 



^(r) = / G{w)e^^^'"-''dudv (1) 

J ~oo 

/oo 
g{v)e-^^-'^-^dxdy, (2) 
-oo 

with w = u&,j. + vBy the spatial frequency vector. The digitization of a signal will invariably 
introduce trains of Dirac distributions. For example, a two-dimensional Dirac train of periods 
li and I2 along e^; and e^, respectively, is defined such that 

00 00 

5{v-v,k)= Yl Hx-ih)S{y-kh), (3) 

i,k=— 00 i,fc=— 00 

with the following Fourier transform relation (see Appendix) 

00 ^00 



i,k=—OD m,n=— 00 



where 



Tik = ihe^ + kkey (5) 

m n 
— e^ H — 
h h 



Hi lb , , 



Another useful distribution is the flat-top window of length A/i (in the e^, direction in 
this case), which we denote by 



-'^In this paper we use a lower case and the corresponding capital letter for a function and its Fourier 
transform, respectively. 
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and the corresponding Fourier transform pair 



, X \ A , . / A , X A , sm (vrnA/i , 

rect -r^ A/i sinc(7rMZ\/J = A/i ^— - — - 

All J nuAli 



3. The Evenly Sampled Grid 

The detection of a signal s (r) from an astronomical source is inevitably achieved through 
a series of transformations. Mathematically speaking, the signal is first convolved with the 
telescope transfer function b (r) such that 



s' (r) = s{r)0b (r) , (9) 

where "(S>" stands for a convolution, while the measured signal t' (r) is a sampled, pixel- 
integrated version of s' (r). For an ESG, the sampling is done in an even manner with a 
Dirac train as defined in equations ([3]) and (jl]). More precisely, for rectangular pixels of 
widths All and AI2 we write 



t'{r) = t{r). J2 Hr-r.,) 



i,k=—oo 

00 

t (r,fe) 5 (r - r,fe) , (10) 



i,k=— 00 



with 



t{r) = s'(r)(g)p(r) 

= [6(r) ®p(r)] ®s(r) (11) 



The convolution 
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h{r) = [b{r)0p{r)] (13) 

stands for what is commonly described as the point spread function (PSF). Using equations 
([2]) and (jlj), and the properties of the Fourier transform for products and convolutions of 
functions, we find that 



T'(w) = T(W)®-^ V 6{W-Wmn) 

m,n=—oo 

CO 

— T{w-w^n), (14) 



m,n=—oo 



with 



T(w) = H{w)S{w 



= 5 (w) P (w) 5 (w) , 



(15) 
(16) 



and 



P (w) = AI1AI2 smc{7iuAl^)smc{7ivAl2) . (17) 

Equation (ITO!) is only valid for the idealized case of an infinite array. In reality this rela- 
tion should be multiplied by an aperture function of appropriate width and shape. Although 
we will take this restriction into account when analyzing the effect of missing pixels in §5.2[ 
we will for the moment simplify our analysis by assuming that the array is sufficiently large 
so that equations ffTOj) and f[T^ are suitable approximations. 



3.1. Interpolation 



The ESG studied in the previous section is the simplest representation that can be 
given for a sampled set of data. As we will see in later sections, we will always seek to 
transform more complicated forms of data grids (i.e., not evenly sampled ones) into ESGs to 
facilitate analysis; this will invariably require the interpolation of sampled quantities from 
different locations. Also, one might inquire about quantities at positions where there are no 
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samples. For example, questions such as "What is the intensity at position A, where there is 
no sample, and how does it compare to the flux at position B, C, and D on this map?" are 
common when analyzing astronomical images. It is, therefore, often necessary to generate a 
new interpolated map from the data set expressed through equation flTU]) . 

Given a weighting function w (r), any value Zi^t (r) to be assigned to an interpolated 
point can be expressed as 



;i8) 



i=l 



where z (rj) is the value associated with the zth of the n data points used for the interpolation. 
The quantity 



n[Y) 



E 

i=l 



W [Y — Yj 



-1 



(19) 



is the normalization factor, which is a function of the position of interpolation. The genera- 
tion of an interpolated map is equivalent to the convolution of the initial data set with the 
weighting function followed by the normalization and re-sampling of the data. This can be 
ascertained through a comparison of equation (ITSl) with 



tint (r) 



5{Y-Y,t - apg) ■ {n (r) ■ [t' (r) ® w (r)]} 

oo 

t (r) ■ (5 (r — r 

i,fc=— oo 

I 

r) ■ ^ t (rjfc) w{y- Yik) 



'■,t= — OD 
OO 

^ 5 (r - r,t - apg) ■ [ n (r 

:,t= — oo 
OO 

:,t= — oo 



i,k=—oo 



where r^^ is deflned as in equation ([5]) and 



h . h 



P 



w [r 



(20) 



(21) 



is the displacement vector specifying the position of the interpolated grid tint (r) in the 
relation to the initial grid. It is important to realize that because of the evenness in the 
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sampling distribution of the original map t' (r) the normahzation factor n (r) will be periodic 
in character with the same periods (i.e., li and I2) as the original sampling Dirac trairj^. As a 
consequence, it will take a common value for all interpolated points similarly located within 
a one-period segment anywhere on the grid (see Appendix). More precisely, data resulting 
from interpolations at points at r and r + rj^, for any integer i and k when is defined as 
in equation ([5]), will have the same normalization factor. Therefore, when the re-sampling is 
done using the same spatial sampling rate as for the original grid (as is the case in eq. [20] ) 
we can write Fourier transform of tint (i") as 



■ 7T T (W - Wmn) , 

m,n=—oo 

(22) 

where c is the constant value associated with n (r) for this particular re-sampling process. 
Equation fl22l) contains multiple copies of T(w), one for each pair of m and n. If the 
Nyquist sampling criterion is satisfied (see Appendix), then the high frequency copies may 
be removed with negligible aliasing by choosing the weighting function such that W (w) ~ 
when |w| > |wmn| /2 (when m 7^ or n 7^ 0). Equation (122|) then simplifies to 



- int 



W 



^ 00 



w 



Wsi)e 



-j27rwsfap 



s,t= 



00 

Tint(w) = — ^ J2 W^(w-w,,)r(w-w,t)e--''2-™-, (23) 

(^1^2) s,t=-oo 

and 



tint (r) = [t (r) ® w (r)] ■ — > 5 (r - r^fc - a^^) . (24) 

nh . , — 

i,k=—oo 

The only difference between tint (j) and t' (r) (see eq. [ID]), besides the overall scaling 
factor and translation, is the presence of the convolution by w (r) for the former. It is 
therefore apparent that w (r) can serve not only as a weighting function for interpolation 
but also as a smoothing kernel, as its effect is functionally similar to that of the PSF h (r) or 
any other function that can be applied to t (r) (see eq. [11]) before or during the sampling 



^It should be noted that the normahzation function wih be constant for an infinite grid when W (w) = 
for |u| > {2h)-' or \v\ > {2h)~' (e.g., sine weighting functions of corresponding widths in normal space) . 
This condition must be strictly enforced in order to obtain a constant normalization factor while satisfying 
the Nyquist sampling criterion. 
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process leading to t' (r). It therefore follows that the weighting functions also possesses 
spectral filtering qualities, as the base spectrum T (w) is multiplied by its Fourier transform 
W (w) (see eq. [23]) . One can, in fact, take advantage of this property in some cases. For 
example, the extraction of a signal from noise can be optimized by matching the spectral 
shape of the Fourier transform of the weighting function (more appropriately named the 
"filter" in this case) to that of the signal itself (if such information is available a priori). 



This is a result commonly established through the so-called matched filter theorem (IHaykin 



19831 ). It is to be noted, however, that optimization of the signal-to-noise ratio through 
filtering is not our goal. As will be made evident in §3.21 besides its fundamental role in 
the interpolation process we are also concerned with determining the effects of the weighting 
function on the spatial resolution of a map. In general, the spatial extent of the smoothing 
kernel will always be significantly smaller than that of the optimized matched filter. 

We now investigate the case of a map resulting from a re-sampling process where we 
seek to increase the density of samples. For example, a map with half the sampling periods 
as the original (i.e., of periods /i/2 and ^2/2) will consist of the combination of four different 
re-sampled maps ti (r), t2 (r), ^3 (r), and (r) that all share the same sampling periods as 
the original map, but translated relative to each other. More precisely, we define 



tl (r) — tint (r) |p,g-+oo (25) 

h (r) = tint (r) |p=2,g^oo (26) 

^3 (r) = tint (r) |p^oo,g=2 (27) 

h (r) = tint (r) |p=2,g=2 • (28) 

In other words, ti (r) is re-sampled at the same positions as the original map while t2 (r), 
ta (r), and t4 (r) are relatively shifted by ^e^,, ye^, and ^e^, + ^e^^, respectively. Calculating 
and summing the corresponding Fourier transforms (using eq. [23]) we find for the combined 
map that 

Ts (w) = Ti (w) + T2 (w) + T3 (w) + T4 (w) 

00 

= 7-7-2 Yl [ci + i-irc2 + i-iyc3 + i-iy^'c,]w{w-wst)T{w-wm) 



where Cj is the normalization constant associated with tj (r) . Correspondingly, we further 
define Acj = Cj — ci, for j = 2, 3, 4, and we rewrite equation ( !29l) as 



- 9 - 



T,(w) = -—2 Yl {4ciiy (w - 2w,t) T (w - 2w,j) 

+ [(-1)^ Ac2 + (-1)* Ac3 + (-1)'+* AC4] Wiw- Wst) T (w - w,t)} . (30) 

As will be soon discussed in §3.2[ the magnitude of the coefficient Acj depends on the 
width of the weighting function w (r); the wider the function, the smaller the coefficient, and 
vice-versa. As we will see below, there are good reasons to limit the width of the weighting 
function, but if we assume for the moment that w (r) is such that Acj ~ 0, then 

T.(w)^— E W^(w-2w,,)r(w-2w,,), (31) 

and 

00 

(r) ^ [t (r) ® «; (r)] . 0- ^ " t) " (^2) 

i,k=—oo 

It is instructive to compare this result with the corresponding equation for an ESG 
t" (r) , similar to that of equation (ITU]) but with half the sampling interval in each direction 



t"(r) = t(r). J2 '^(^-t) ^^^^ 

j,fc=— 00 

4 °° 

^"(^) = TT ^ T(w-2w.,). (34) 

^ ^ s,t=-oo 

Again we see that, apart from a multiplication factor, the (approximate) interpolated 
grid tg (r) differs from t" (r) by the presence of the convolution by w (r). Although equations 
( I3T]) and f l32i) are approximations and the aforementioned correspondence between ts (r) and 
t" (r) may fail for a given weighting function (i.e., Acj 7^ in general in eq. [30]), this 
simplification is often reasonable (see §3.21) . At any rate, one can more closely approach the 
idealization of equation (!32|) by broadening the width of the weighting function w (r) (with 
a corresponding loss in spatial resolution, however). 
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3.2. Selection of the Weighting Function 

There is a fair amount of subjectiveness in choosing the specific form and characteristics 
of a weighting function. We choose the following two criteria as guidelines for achieving this: 

1. The function must be sufficiently broad so that its amplitude is large enough at the 
interpolated positions, while not being too broad to significantly degrade the resolution 
of the map. 

2. Its spectral extent must be such that it filters out spatial frequencies for which \u\ < 
(2/i)~^ or \v\ < (2/2)"^ (see eq. [22] and the discussion that follows). 

We now show how this can be practically implemented by considering the case of SHARC-11 
(or SHARP) where Ali = li = AI2 = h- Furthermore, we approximate the SHARC-II PSF 
with the following Gaussian profile 



h{r) 
H(w) 



1 I'M 

2 \ cr 



■ |w 



The PSF size is usually defined by its full-width -half- magnitude (F WHM) , which is 



approximately 9 arcseconds for SHARC-II at 350 /xm (IDowell et al.l 120031 ). This gives a 



FWMH/a/8 In (2) ~ 3.8 arcseconds; we will use standard deviations to specify widths of 
Gaussian PSFs. With this definition, the one-sided bandwidth associated with SHARC-II 
at 350 fim is 



2na 23.9 

Since li = I2 — 4.7 arcseconds, then > 2/i and the Nyquist sampling criterion is 
met, as previously assumed. If we were to choose the weighting function w (r) to also be 
Gaussian, and of width zo, then to satisfy Criterion 2 above we must have 

zu > — ~ 1.5 arcseconds. 



Taking the lower limit for the size of the kernel, we can evaluate the new resolution of 
the map with 
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/ p 

Vo"^ + tA7^ = a\ 1 -\ — ^ 1.07(j = 4.1 arcseconds, (35) 
V n'^a'^ 

which corresponds to an equivalent PSF width of 9.6 arcseconds for SHARC-II and a loss 
of approximately 7% in spatial resolution. Finally, the relative amplitude of the weighting 
function at a distance of one-half pixel away from the position of interpolation would be 



27rzu'^w (r) 



— ii 

^~ 2 



0.29. (36) 



Although equations and satisfy Criterion 1 above and one could reasonably 
choose the corresponding weighting function to interpolate a map, one should nonetheless 
verify that the coefficients Acj resulting from the interpolation process (see §3.11) are suffi- 
ciently small when seeking to increase the density of samples in the final grid. Doing so will 
ensure that the approximation leading to equation (1521) is adequate, for example. Figure [T] 
shows a map (top) of the normalization function n (r) for a SHARP ESG with the lower-limit 
weighting function considered above {tu = li/n ~ 1.5 arcseconds). The top most curve in 
the lower part of the figure is a cut through a row or column of pixels for the normalization 
map. The bottom two curves are similar cuts for weighting functions of = 1.3/i/7r and 
l\l \fv: ^ 1.8 /i/vr, respectively. It is clear from these curves that the relative amplitude of the 
Acj coefficients, which can be asserted from the level of ripple on the curves, exhibits a strong 
dependency on the width of the weighting function. For example, the two larger weighting 
functions in the lower part of Figure [T] exhibit variations in amplitude of 11% and 1% for a 
loss in spatial resolution of 13% and 22%, respectively. This behavior is traced to the fact 
that a larger weighting function will more completely cover the space located between neigh- 
boring sampling positions, hence the existence of a smoother normalization function n(r). 
This will be better visualized for the general case with the graph shown in Figure [2] where 
trends in normalization function (solid line) and spatial resolution degradation (dashed lines) 
with smoothing kernel size are plotted (see the corresponding caption). 



4. The Regularly Sampled Grid 

We define a Regularly Sampled Grid (RSG) as being a generalization of the ESG dis- 
cussed in the previous section. That is, a RSG has a well defined periodicity (just as the 
ESG), but it is a grid for which the pattern of Dirac distributions is more complex. While 
along a coordinate axis of the ESG there is only one Dirac distribution for a given period, a 
RSG may have many Dirac distributions (not necessarily evenly spaced) over the same inter- 
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val. This difference is illustrated in cases (a) and (b) of Figure [3] for one-dimensional versions 
of an ESG and a RSG, respectively. Practically speaking, a RSG would be encountered any 
time that maps of similar characteristics, but translated relative to one another, are com- 
bined together to form a unique and final map. An example for this would be astronomical 
images of a given object at different pointing positions. 

It should be apparent from this discussion and Figure ^ that a RSG, which we again 
denote by t' (r), can be simply expressed as a combination of a set of relatively displaced 
ESGs with 



^ (5 (r - Tifc - dp) 



(37) 



p=i 



_i,fc=— oo 



where t (r) is defined in equation ffTTl) and 



dp (r) = XpS^ + ypSy (38) 

is the relative displacement associated with the pth of the Ug ESGs that make up the RSG. 
It is straightforward to calculate the Fourier transform of equation (I57|) to get 

rig oo 

T'(w) = — V V r(w-w„„)e-^2— ■'i^ (39) 

p=l m,n=—oo 

with defined in equation 

The important aspect to emphasize for the interpolation of a RSG is that, as was the case for 
an ESG, the normalization factor n (r) in equation (ITSl) is common to all interpolated points 
similarly located within a one-period segment anywhere on the grid. This is illustrated with 
the vertical broken lines in Figure [31 The existence of such a common normalization factor 
could be effectively adopted as the definition for a RSG. 

If we interpolate our RSG using a weighting function w (r) that satisfies Criterion 2 
above, then the high spatial frequency components of the spectrum will be filtered out. 
Therefore, starting from equation (jSHD, and using steps similar to those that led from equation 
([H]) to equations and fl25]) . the resulting interpolated grid becomes 

oo 

^int(w) = -^ J2 W^(w-w,,)r(w-w,,)e-^'2----^-, (40) 

s,t=-OC 
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and 

oo 

tint (r) = [t (r) ® w (r)] ■ ^ V 6 {r - nk - a^g) . (41) 

i,k=—oo 

Once again c and a.pg denote, respectively, the common normalization factor and the 
displacement of the new interpolated grid in relation to the original grid (see eq. |21j). 
Obviously, the same comments apply for the map resulting from the re-sampling of a RSG 
here as for an ESG in ^ That is, the most notable effect of the interpolation/re-sampling 
process is the presence of the convolution by the weighting function in equation (HTj) . 

5. The Irregularly Sampled Grid 

An Irregularly Sampled Grid (ISG) can manifest itself in different ways. For example. 
Figure [3]c shows a case where the distribution of Dirac functions within a given base period 
(of length / in the figure) is not the same from one interval to the next. Another possibility 
is shown in Figure [3]d where no Dirac distributions are present for some intervals. This can 
be likened to situations where pixels are missing from an array detector (see below). 

The problem in the analysis of an ISG is twofold. First, there is no simple way of 
expressing the Fourier transform of irregularly spaced Dirac distributions such that the 
spectrum will show a repeating pattern of some frequency, as is the case for an ESG or 
a RSG. Moreover, the lack of regularity in the positions of the Dirac distributions implies 
that there does not exist a common normalization factor when performing interpolations to 
create an ESG from an ISG (see below). Nevertheless, it is still possible to analyze some 
specific types of ISGs. We deal with two possible cases in what follows. 

5.1. The Combination of Relatively Translated and Rotated ESGs 

It often happens that an astronomical source will be observed at different times, when 
it is at different locations and orientations on the celestial sphere. Invariably, we seek to 
combine the resulting images to form a final map of the object. If the array detector (which 
we assume perfect and therefore able to generate ESGs of data) used to record the images 
is part of an instrument that is unable to precisely track the apparent rotation of the source 
on the sky, then the different images of the source will be sampled with ESGs that will be 
rotated and possibly translated relative to each other. A simple example is shown in Figure 
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m where two ESGs are combined: one rotated by 10 degrees with respect to the other. These 
grids are not relatively translated (see the caption). 

Perhaps the most important aspect of Figure H] is the fact that the combination of the 
two ESGs produces a grid which has an irregular pattern of Dirac distributions, as can be 
asserted by the coverage of a predetermined weighting function (shown with large empty 
circles in the figure). Clearly, any weighting function w {r) is likely to cover a different 
number of samples at different positions on the map. The main consequence resulting from 
this fact will be the absence of a common normalizing factor at the different locations where 
interpolations are performed (note that the normalization function in eq. [19] will not be 
periodic). We can therefore expect that interpolated maps originating from JSCs will be 
more complex than those resulting from ESGs and RSGs. 

Another point to consider is the possible relative rotation between the different maps to 
be combined. Because of this, we will do well to use the fact that the Fourier transform of 
a rotated map is the rotated Fourier transform of the original (i.e., without rotation) map. 
That is, if we have the Fourier pair 



where R stands for the rotation operation (i.e., matrix). Because of this property of the 
Fourier transform we can express an ISG t' (r) composed of Ug rotated and translated ESGs 
and its Fourier transform as 



g{r)^G (w) , 



then it is also true that (see Appendix) 



g (Rr) ^ G (Rw) , 



n, 



■9 



oo 




(42) 




T(w 




(43) 



where Hp and dp = Xpe^ + Upey are, respectively, the rotation matrix and the translation 
vector corresponding to grid p. Just as for the RSG our goal is to generate an ESG tint (r) 
of periods li and I2 (along the x and y axes, respectively) from t'(r). Although, as was 
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previously pointed out, we cannot express the interpolation process with a simple convolution 
with a weighting function w{r), we can still use the general expression given in equation 
fl2U]) . That is, with denoting the origin of the new interpolated grid (see eq. [2T]) we 
have 



^int (r) 



[t' (r) (g) w (r)] • n (r) ^ 5 {r - r^fc - a 

l,fc= — oo 

rig oo 

i (r) • ^ ^ 5 (r - RpFifc - 



p=l i,fc=— oo 



to r 



■n r 



5 (r - Fifc - 



i,fc=— oo 



Calculating the Fourier transform of equation (jH]) we get 



(44) 



Tint (w) 



^ rig oo 

TeJ^YI T{w-RpWst)e 

p=l s,t=—oo 
^ oo 



W(w) > ® iV(w) 



m,n=—oo 



which, using the assumption that W (w) is such that it filters out the higher frequency 
components of the spectrum, can be approximated to 



Tint (w) 



^ e-^'^'^^™-^^^ {[T(w) Vr(w)] ® A^(w)} 



w=w— w„ 



m,n=— oo 



Correspondingly, we can approximate equation (l44l) to 



(45) 



tint(r) = {[t(r) ®w(r)] ■n(r)} ■ ^ V 5 (r - r^fc - a. 



i,k=—oo 



(46) 



Equation fl45|) shows best the effect of the lack of a common normalization factor on 
the interpolation process. Since the Fourier transform N (w) of the normalization function 
is convolved with the weighted (and low-pass filtered) spectrum T (w) W (w), the resulting 
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spectrum of the interpolated ESG tint (r) is broadened by A^(w). It is interesting to note 
that w (r) and n (r) have opposite effects on the signaL That is, the weighting function 
restricts the extent of the spectrum, while the normalization function extends it. 

Given such an ISO, it should be in principle possible to evaluate n (r) and quantify 
its effect. In particular, one should ensure that the two previous criteria (see §3.2p used to 
select the weighting function are met. Optimally, n (r) will be sufficiently slowly varying, 
and of low enough amplitude, that the spectral broadening will be minimal. To make this 
clearer, we show in Figure an example consisting of a combination of two relatively rotated 
ESGs, similar to those of Figure H] (i.e., one rotated by 10 degrees with respect to the other, 
and no relative translation between the two). The map at the top of the figure is for the 
normalization function n (r) of the resulting SHARP ISG using a weighting function with 
w = Ii/ti {li = I2 — 4.7 arcseconds for SHARP). The top most curve in the lower part of the 
figure is an arbitrary cut through a row of pixels for this normalization map. The bottom 
two curves are similar cuts for weighting functions of w = l.S/i/vr and l\l ~ I.S/i/tt, 
respectively. The black dots highlight the values taken by n (r) for a re-sampling onto an 
ESG at the original sampling rate. It is clear from the top two cuts that the normalization 
factor is not constant in general. One can also assert from this that the spectrum due to a 
broad source (in relation to the size of the map) would be significantly more broadened by 
the weighting function that produced the top curve (i.e., with w = h/n) than by the other 
two. 



5.2. The Effects of Missing Samples 

It is a common, if unfortunate, fact that detectors arrays used in astronomy will often 
contain pixels that are either performing significantly below specifications or are completely 
unusable. Astronomers usually work around the difficulties occasioned by these so-called 
"missing" pixels by dithering the array during observations, thus ensuring complete mapping 
of the source under study. It would be, however, instructive to analyze and quantify the 
impact that missing pixels would have on the representation of astronomical signals without 
such corrective techniques. 

Although we will take into account the fact that the map obtained from the array 
is composed of a finite number of samples, our approach will consist of first temporarily 
lending it an infinite character and then removing it. More precisely, although the size of the 
(rectangular) detector array considered in this section is Nili x N2I2, we first assume that 
the two-dimensional pattern of N1N2 pixels (including the missing pixels) repeats infinitely 
in all directions. The underlying assumption is that the finiteness of the map will be restored 
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in the end by windowing with the appropriate aperture function. Using this approach we 
express the sampled signal (before windowing) as 

oo 

t'(r)=t(r)-5^ 5 (r - d,, - rpij , (47) 

pix s,t=—oo 

with Tpix the position of a pixel on the array and (take note of the periods) 

dst = sNihe, + tN^hey. (48) 

That is, we first account for the pixels of the finite array through the summation Xlpix' 
and then associate a Dirac train of periods {Nili, N2I2) to each pixel. These Dirac trains 
are relatively translated in space and are accounted for by the summations on s and t in 
equation fHTl) . The Fourier transform of the combined Dirac trains is 



5{r-dst- Tpix) ^ 



pix s,t=—oo 



pix 



-j27rw-rpix 



^ 00 



(49) 



m,n.=— 00 



with 



m n 

^rnn = ^e. + ^e,. (50) 

Note that the minimum separation between two Dirac distributions in frequency space 
is 1/Nl, where Nl is the greater of Nili and N2l2- We write the right hand side of equation 
din]), which we denote by D (w), as follows 

^ (w) = — y2 E (W„„) 5 (W - ^rnn) , (51) 
m,n=—oc 

with 

^ ^'^'""^ = S e-^-2™™-p-. (52) 

pix 

The effect of missing pixels can now easily be quantified, at least in principle, by omitting 
them from the ^pj^ summation in equations fl52l) . As an example, consider the case where 
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the M X M pixels in the "top-right" comer of the array are missing. For this particular 
example, we will do well to make the following substitutions 



pix i k 

where i and k indices stand for the columns and rows of the detector, respectively. We can 
then transform equation fl52l) to 



Afi-l 7V2-Af-1 

j2nkn/N2 



e 

i=0 k=0 
Ni-M-1 N2-I 
_^ ^ ^-j2TTim/Ni _ ^ ^-j2nkn/N2 

i=0 k=N2-M 
j-Kn{M+l)/N2 



(_!)"+" c'^^'^IN, sin (7rm) sin [vm {N2 - M) j 



N1N2 [ sin(7rm/A^i) sin(7m/iV2) 

I / .xm ,-^^(M+i)/jVi Sin[7rm(A^i-M)/iVi]sin(7rnM/Ar2) \ 
^ ^ sin(7rm/Ari) sin(7™/Ar2) J" ^ ^ 

When M = (i.e., when all the pixels are accounted for) equation fl53|) simplifies to 



1, m = mWi andn = n'N2 
0, elsewhere 



with m' and ra' some integer numbers. That is to say, equation (1511) then becomes 



00 



m,ri=— 00 

which is, as it should be, the same result as was obtained earlier with equation (jl]) for the 
ESG. 

Although it is necessary to plot E [wmn) to assess the effect of an arbitrary distribution 
of missing pixels, it should now be clear from equation fl53l) that its amplitude is in general 
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non-zero for all values of m and n. In other words, by removing even only one pixel we 
went from a case where we only had Dirac functions at frequency intervals of and 
to a situation where they are separated by intervals of only (A^i/i)~^ and (A^2^2)~^- It is 
important, however, to quantify the relative magnitude of these Dirac distributions. For 
example, returning to equation (!53l) pertaining to our case of the M x M missing "top- 
right" pixels, we find that 

E (wio) 
E (woo) 



This last relation yields the perhaps intuitive result that when only a small number 
of pixels are missing the amount of contamination determined by the ratio expressed in 
equation fl5^ is approximately equal to the ratio of the number of missing pixels to that 
of good pixels. We should, however, resist the temptation to generalize this result, since 
different distributions of missing pixels would give different levels of contamination. 
Especially if they are not concentrated in one part of the array, as is the case here. An 
example is shown in Figure [6] where the function E (wm^„) is plotted for three different cases. 
Starting with the nominal SHARP 12 x 12 array, E (wmn) is shown for n = when no 
pixels (black curve and dots), ii) 16 randomly positioned pixels (red curve and dots), and 
in) the 4x4 "top-right" corner pixels (blue curve and dots) are missing. Contrary to the 
case of an ESG (corresponding to the black dots) where E (wm„) = when \m\ ^ 0, 12, 
an ISG (red and blue dots) will in general have E (wmn) 7^ for all m and n. It should 
be clear that E (w^n) acts as a mask that will or will not allow the appearance of Dirac 
distributions that are more closely spaced in frequency depending on whether or not there 
are missing pixels in the array. This serves to emphasize the fact that missing pixels will 
bring some spectral contamination (i.e., aliasing) in the sampled signal. 

This becomes more evident if we calculate the spectrum of the measured signal from 
equations (H7j) and (ISTIl 

^ 00 

T' {^) = Yj- E (wm„) T (w - w™„) . (55) 

m,n=— 00 

We see that the different rephcas of T (w) are spaced in frequency according to equation 
f l50|) (compare this with the case of the ESG in equation ([6]) where the spacing between 
replicas is A^-times greater). Contrary to the case of an ESG, a convolution with the usual 



M I sin [tt [Ni - M) /Ni] / sin (vr/A^i) | 



N1N2 - M2 
forM < A^i. 



(54) 
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weighting function w (r) will not restore a low-pass filtered version of the t (r) map. To make 
this clear, we first define a new function Y (w) such that 

^ (w) = ^ (wm„) T (w - w^„) . (56) 

Then proceeding with the usual interpolation defined in equation (ITSj) we find that 

oo 

7^int(w) = -— ^ E e-''''^"^*-^-{[r(w)iy(w)]®iV(w)}__^^ (57) 

and 

^ oo 

tint(r) = {[y(r)®^i;(r)]-n(r)}-— V 5 (r - r^fe - ap^) , (58) 

2,fe = — oo 

where S-pq is defined in equation (I2T1) and denotes, once again, the position of the origin of 
the interpolated grid. These equations show that maps resulting from the interpolation of 
ISGs containing missing pixels suffer from both spectral aliasing (from the presence of Y (w) 
in lieu of T (w) in eq. [57]) and broadening (because of the presence of (w)). 

We once again stress the realization that missing pixels will bring some amount of 
aliasing that will he impossible to remove with a reasonably sized weighting function. The 
concept of Nyquist sampling can even lose much of its meaning and usefulness in a situation 
where too many pixels are missing, or when the level of contamination due to spectral aliasing 
is comparable to the noise level present in the map. This fact strongly underlines the necessity 
of performing adequate dithers or other scanning strategies when observing with an imperfect 
detector array. 

We complete the analysis by performing the windowing mentioned earlier, which trans- 
forms equation flS^ to 



t'int (r) 



^int (W) 



''int 



rect 



X 



rect 



y 



Tint (w) ® [A^i/iiV2/2 sine (TruNih) sine (vrf A/'a/a)] • 



(59) 
(60) 



For reasonably large detector arrays we do not expect that the presence of the sine 
functions will be of any significance due to their spectral narrowness relative to the extent 
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of T (w) (or Y (w) W (w)). Because of this, our periodic depiction of the array, upon which 
our analysis rests, is justified. 

6. Summary 

In this paper we addressed the question of astronomical image processing from data 
obtained with array detectors. We defined and analyzed the cases of evenly (ESG), regularly 
(RSG), and irregularly (ISG) sampled grids for idealized and realistic detectors. We focused 
on the effect of interpolation on the maps, while using a Gaussian kernel to accomplish this 
task. In all cases (i.e., ESG, RSG, and ISG) we have applied the method of weighted averages 
(eq. [H]) to produce a map interpolated on a finely spaced grid. 

We defined an ESG as a map where the signal to be analysed is digitized with a simple, 
two-dimensional, train of Dirac distributions evenly separated with well-defined spacings 
(see eq. [3]). Moreover, since the ESG is the simplest way to represent and analyse a set 
of sampled data, we always sought to transform a non-evenly sampled grid (i.e., a RSG 
or an ISG) to an ESG through the process of interpolation. While studying the ESG, we 
found that the interpolation process invariably leads to a loss in spatial resolution; this loss 
grows with increasing width of the smoothing kernel (this result is true in general, i.e., when 
considering RSGs and ISGs). When an ESG is re-sampled at the same rate as the original 
map the interpolation process can usually be adequately taken into account by replacing 
the original signal by its convolution with the weighting function (see eq. [21]). The same 
is not true in general, however, when the final ESG is re-sampled at a rate different than 
that of the original map (see eq. [SO])- This is due to the fact that the normalization 
function that is intrinsic (and necessary) to the interpolation process is not constant but 
a function of position (although it is periodic with periods corresponding to the sampling 
rates). The aforementioned replacement of the original map in the interpolation process by 
its convolution with the weighting function will only be adequate in such cases when the 
latter is sufficiently broad relative to the spacing between samples (see Figs. [Hand [2D. 

We defined an RSG as a generalization of an ESG such that it consists of a combination 
of a number of relatively translated ESGs of similar sampling rates. All the results obtained 
for the ESG can be generalized to the RSG. 

We analyzed two different types of ISGs: the combination of relatively translated and 
rotated ESGs, and ESG-like maps with missing samples (e.g., data grids made with detectors 
exhibiting dead pixels). In the first case, the interpolation process cannot be simply repre- 
sented by a convolution of the original signal with the weighting function, as this operation 



22 



must be subsequently multiplied by the normalization function (see eq. [16]). Because of the 
irregular nature of the new map, the normalization function is not periodic in general and 
will add structure to the spectrum (i.e., the Fourier transform) of the map. More precisely, 
the spectrum of the source (filtered by the spectral profile of the weighting function) will be 
broadened through its convolution with the Fourier transform of the normalization function. 
This effect is a function of the size of the weighting function, as the spatial variation of the 
normalization function grows larger for smaller kernel widths. Although these results also 
apply to maps exhibiting missing samples, we found that these further suffer from spectral 
aliasing that may reduce or negate the usefulness of the Nyquist sampling criterion in ex- 
treme cases. This fact strongly underlines the necessity of performing adequate dithers or 
other scanning strategies when observing with an imperfect detector array. 
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In this Appendix we provide a few simple derivations to justify some of the results used 
in the text. 



The Fourier transform of a Dirac train can easily be determined by first calculating the 
associated Fourier series. In the one- dimensional case we have 



A. Appendix 



A.l. Fourier Transform of a Dirac Train 




(Al) 



i=—oo 



n=—oo 



since the Fourier series for a periodic function g (x) of period / is defined as 
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n=— oo 

with the Fourier coefficient G [n] 

1 M 

G{n) = - / g{x)e-^^^''Tdx. 

Before calculating the Fourier transform of equation flAip . we consider the so-called 
duality property of the Fourier pair in equations ([T]) and ([2]). More precisely, we mean that 
if for a function / (r) we have 

/(r)^F(w), 
then it must also be true that for F (— r) we have 

F(-r)^/(w), 

as can be readily verified by inspection of equations ([T]) and ([2]). It follows from this that 
since 

5(r-ro) 

then 

gi2uwo r ^S{w-Wo). (A2) 

Using the one-dimensional version of equation flA2l) to calculate the Fourier transform 
of equation flAip we find that 

oo ^ oo 

J2 s{x-zi)^- J2 ^(«-y)- (A3) 

i=— oo ri=— oo 

The two-dimensional generalization of this result is straightforward and leads to equa- 
tions da]) and (gD . 



A. 2. The Nyquist Sampling Criterion 



The Nyquist Sampling Criterion can be understood with equation (IA3P and the prod- 
uct/convolution property of the Fourier transform. This property states that the following 
Fourier pair is valid 

f(r)g(j)^F{w)®G{w) 

for two functions / (r) and g (r) , as can easily be verified from the definition of the Fourier 
transform. Therefore, for the sampling of a function t [x) we have 



oo _ oo 

t{x)- ^ 5{X-Il) ^ T(u)®y ^ ^^M-y) 



^ oo 



n 



which implies that the base spectrum T {u) is repeated in frequency space at an interval equal 
to the sampling period of l^^. It is apparent that in order to avoid any cross-contamination 
between the different spectral replicas of T (u) the following relation must be enforced 



T{u)=0 for \u\ > — (A4) 

Relation flA4p is the one-dimensional mathematical equivalent of the Nyquist sampling 
criterion, which states that in order to recover the base spectrum T [u] from its sampled 
version (through spectral filtering) the sampling frequency must be at least twice as large as 
the frequency extent of T {u) . 



A. 3. Normalization Factor 

We know from our analysis that an interpolated map tint (r) resulting from a previously 
sampled data set t' (r) is given by 

oo 

itint(r)= (5(r-r,t-apg)-{n(r)-[t'(r)®M;(r)]}, 
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which we transform shghtly to 
^int (r) = 

As usual n (r) and w (r) are the normahzation and weighting functions, respectively, and 
the vectors Vst and S-pq are given by equations ([5]) and (12T1) . When the original map is an ESG 
the normalization function is periodic and can therefore be expanded with a two-dimensional 
Fourier series 

oo 

n (r) = ^ N (z, k) (A6) 

i,fc=— oo 

with N {i, k) the corresponding Fourier coefficient and Wj^ given by equation (j6]). 

From equations (]A2p and (]A6P we can write the Fourier transform of equation (!24l) as 



n(r] 



5 fr - r^j - a, 



[t' (r) ® It; (r)] 



(A5) 



s,t= 



but since 



^ N{t,k)S{w-w^k) 

i,k=~OD 

® [T' (w) W (w)] 

^ oo oo 



^ oo 



m,n=— oo 



w - w„„ - Wife) e 



W - Wm.r,. e 



-j27r(w-Wife)-ap 



-j27rw-ap 



i,fc=— oo 



m,n=—oo 



® [T' (w) (w)] 

oo oo 

— J2 AT (i. A;) e-^'2""''^ "^« J] 5 



i,fc=— oo 



m,n=—oo 



[T (w)W^(w)], 



^ 5 (W - Wrnn - ^ik) = - ^m'n') 

m,n=—oo m'n'=—OQ 



with m! = m + i and = n + then 
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Tint fw) 



J2 N{z,k)e-^ 

,i,k=—oci 

[T' (w) W (w)] 



27rw,-t - a- 



ik '^pq 



hi 



— ^ 5 (w - Wm'n') e'^ 



27rw-ai, 



m' ,n '= — oo 



However, we can use equation flA6l) one more time to transform this last relation to 



Tint fw) = n fa, 



■PQJ 



^ oo 



w - w™,„ e 



-j27rw-ap 



m,n=— oo 



[r'(w)i^(w)], 



and 



tint (r) = n (apg) ^ 5 (r - r^i - ap^) ■ [t' (r) (g) w (r)] 



(A7) 



s,i=— CO 



Evidently n (apg) is constant for a given interpolated map, but will vary as a function 
of the displacement of the new sampling grid (defined by a.pg) relative to the original one. 
Equation (1A7I) leads to (and justifies) equation (12^ provided that we set c = n (apq). 



A. 4. Fourier Transform of a Rotated Map 

Finally, we prove the result used in ^that the Fourier transform of a rotated map is 
the rotated Fourier transform of the unrotated map. To do so we subject a two-dimensional 
map g (r) to a rotation R and calculate the Fourier transform of the transformed map g (Rr) 
with 



G' (w) 



g (Rr) e-^^^'^-'dxdy 



We now make the change of variable r' = Rr to get 



G' (wl 



/oo 
-oo 

/oo 
g (r') e-'^^^^-^^-^'dx'dy', 
-oo 



(AJ 
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where the last transformation was made possible by the fact that the inverse of a rotation 
matrix equals its transpose. We therefore find from equation flASp that if 

g{r)^G (w) , 

then 

g (Rr) ^ G (Rw) . 
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Fig. 1. — A map (top) of the normalization function n (r) for a SHARP ESG with a weighting 
function with zu = li/n (/i = /2 — 4.7 arcseconds for SHARP). The top most curve in the 
lower part of the figure is a cut through a row or column of pixels for the normalization 
map. The bottom two curves are similar cuts for weighting functions of zu — I.SIi/tt and 
h/V^ — I.SZi/tt, respectively. 
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Fig. 2. — Trends in normalization function (solid line) and spatial resolution degradation 
(dashed lines) with smoothing kernel size. The kernel sizes on the abscissa are given as the 
Gaussian width (zu; lower axis) and FWHM (= a/8 In (2) rxj; upper axis), both in units of 
the array pixel separation. The amplitude of the normalization function's spatial variation 
is given as a percentage of the function's average value (see Fig. [T]). The corresponding loss 
in spatial resolution is described by equation fl35|) and the following text. This is plotted 
here for different beam sizes and indicated in units of pixel separation. For example, the 
case of SHARP, with Beam = 9"/4.7" ~ 2, closely corresponds to the second dashed curve 
(from the top). 
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Fig. 3. — One-dimensional examples of (a) an Evenly Sampled Grid (ESG), (b) a Regularly 
Sampled Grid (RSG), and (c) and (d) two Irregularly Sampled Grids (ISGs) are shown. Inter- 
polations at the positions of the two vertical broken lines would require weighting functions 
that have a common normalization factor for the ESG and RSG, but different normalization 
factors for the ISGs. Dirac distributions are shown as vertical arrows. 
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Fig. 4. — A combination of two relatively rotated ESGs. Every small dark dot corresponds 
to a Dirac distribution, and the position of the small empty circle is the origin of the maps 
and of the rotation for one of the two grids. Its relative rotation is of 10 degrees with respect 
to the coordinate axes (also shown). Neither ESG is translated. The four large empty circles 
correspond to the footprint of a predetermined weighting function. 
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I . . . , . . . , . . . I 

Fig. 5. — A combination of two relatively rotated ESGs, as those of Figure HI We show a 
map (top) of the normalization function n (r) for the resulting SHARP ISG using a weighting 
function with = li/n (/i = ^2 — 4.7 arcseconds for SHARP). The top most curve in the 
lower part of the figure is an arbitrary cut through a row of pixels for this normalization 
map. The bottom two curves are similar cuts for weighting functions of tu = l.S/i/vr and 
h/V^ — l.S/i/vr, respectively. The black dots highlights the values taken by n{r) for a 
re-sampling onto an ESG at the original sampling rate. It is clear from the top two cuts that 
the normalization factor is not constant in general. 
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Fig. 6. — The function E (wmn) for n = when no pixels (black), 16 randomly positioned 
pixels (red), and the 4x4 "top-right" corner pixels (blue) are missing. Contrary to the case 
of an ESG (corresponding to the black dots) where (w^o) = when |m| ^ 0,12,..., an 
ISG (red and blue dots) will in general have E (w^n) 7^ for all m and n. The amplitude of 
El (wmn) for the relevant (i.e., integer) values of m are shown by the colored dots. The curves, 
which include computations at intermediate values of m, are only shown to emphasize the 
fact that E iy^rnn) can be interpreted as a mask function (see text). 



